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The so called quantum spin Hall phase is a topologically non trivial insulating phase that is 
predicted to appear in graphene and graphene-like systems. In this work we address the question 
of whether this topological property persists in multilayered systems. We consider two situations: 
purely multilayer graphene and heterostructures where graphene is encapsulated by trivial insulators 
with a strong spin-orbit coupling. We use a four orbital tight-binding model that includes the full 
atomic spin-orbit coupling and we calculate the Z 2 topological invariant of the bulk states as well 
as the edge states of semi-infinite crystals with armchair termination. For homogeneous multilayers 
we find that even when the spin-orbit interaction opens a gap for all the possible stackings, only 
those with odd number of layers host gapless edge states while those with even number of layers are 
trivial insulators. For the heterostructures where graphene is encapsulated by trivial insulators, it 
turns out that the interlayer coupling is able to induce a topological gap whose size is controlled by 
the spin-orbit coupling of the encapsulating materials, indicating that the quantum spin Hall phase 
can be induced by proximity to trivial insulators. 

PACS numbers: 73.22.Pr, 73.43.Cd 


I. INTRODUCTION 

In their seminal papers,® Kane and Mele established 
the existence of two fundamentally different types of band 
insulators with time reversal symmetry in two dimen¬ 
sions, dubbed as trivial and topological. Remarkably, it 
was predicted that monolayer graphene would be topo¬ 
logical, giving rise to protected chiral gapless edge states. 
Importantly, this opened a new venue in condensed mat¬ 
ter physics, the quest of searching and designing topolog¬ 
ical states in two dimensional systems. 

The nature of the topological state in graphene comes 
from the intrinsic spin-orbit coupling (SOC). In partic¬ 
ular, SOC will open gaps of opposite signs at the two 
Dirac points, in contrast with the trivial gap that a stag¬ 
gered potential opens in the honeycomb lattice, with the 
same sign at the two valleys. This twisting of the wave 
functions in the reciprocal space leads to the appearance 
of in-gap states at the boundaries of the material. Sub¬ 
sequent work 2 !® found that the size of the SOC gap in 
graphene was very small, and the attention shifted to 
other systems, such as CdTe/HgTe quantum wells 5 in 
which the quantum spin Hall (QSH) phase was found 6 
as well as to bulk systems, for which the notion of topo¬ 
logically non-trivial insulators was extended. Experimen¬ 
tal evidence for quantum spin Hall phase has also been 
found in other systems, such as Bi(l 11) atomically thin 
layerd®, and InSb/GaSb quantum well d 9 l 1Q [ 

Multilayers of two dimensional materials are also po¬ 
tential candidates to sustain topological states. In par¬ 
ticular, their appealing comes from the the tunability 
of stacking a different number of layers, or even differ¬ 
ent materials. In the present work we will focus on 
the study of a particular type of multilayer systems, 
whose basic building blocks are graphene-like systems. 
We will study mainly two families of multilayers. First 


we consider multilayered systems formed by graphene¬ 
like insulators using the SOC as a free parameter, so 
the main concepts should be suitable for systems such as 
graphene, Silicene 11 !®, Germanene 1 ® or Stanene 15 ! (in 
fact our methods make it easy to ex tended this kind of 
analysis to the case of BismutlPEE 16 , and metal-organic 
frameworks 17 23 !). Second, stacks formed by a layer of 
graphene encapsulated by some trivial insulator with a 
strong SOC. 

From a practical point of view, several reasons mo¬ 
tivate this work. First, there is a generic interest in 
the possibility of engineering the electronic properties of 
two dimensional crystals, such as graphene, h-BN and 
transition metal dichalcogenides, by combining them into 
multilayeriPM26] slacking monolayers of the same type 
is also a very interesting and widely studied possibility. 

Our second motivation is to study the behavior of the 
topological gap as we increase the number of layers in the 
system. In the case of graphene, it is well known that key 
electronic properties, such as the pattern of Landau levels 
and the density of states at the Dirac point are drastically 
modified for bilayeJ^and trilayer graphene.®®. Recent 
experimental work shows that some sort o f magnetic or¬ 
der can occur, even at B = 0 in bilayei®^®, trilayei^ 5 and 
even tetralayerd^. These last trivial symmetry breaking 
states, will compete with the potential topological states 
studied in our work. 

A third motivation comes from recent experimental 
that report a very large enhancement of the spin Hall 
effect for graphene deposited on top of WS 2 , as trivial 
semiconductor with a quite large SOC. This inspires the 
calculation for graphene placed between two insulators 
with a trivial band gap, large SOC and broken inversion 
symmetry, to mimic the properties of WS 2 and related 
transition metal dichalcogenides. 

Furthermore, we have also a formal motivation. It is 
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not obvious a priori that the original second-neighbor 
hopping HamiltoniarP can be applied to multilayer 
graphene 3 ^. In monolayer graphene the p z orbitals are 
strictly decoupled from the s, p x , p y orbitals, due to mir¬ 
ror symmetry with respect to the plane. In the mono- 
layer, SOC mixes p z with p x and p y orbitals of oppo¬ 
site spin and, when treated perturbatively, leads to an 
effective Hamiltoniarf 1 with a spin-dependent effective 
second-neighbor hopping between p z orbitals that con¬ 
serves S z . In multilayer graphene this is no longer true, 
since electrons in a p z orbital in one layer can hop to 
the s orbital of atoms in the next layer. When SOC is 
added to the model, we expect that this s — p z mixing 
naturally leads to spin-mixing terms in the Hamiltonian, 
which is indeed the caseP^ The presence of this spin- 
flip channel interaction cast a doubt on the validit y of 
the spin-conserving Kane-Mele model for multilayer ^ 40 41 
and motivates our choice of the standarcPE^HUl f our or _ 
bital tight-binding calculations. 

The rest of this work is organized as follows. In sec¬ 
tion |II| we briefly review the tight-binding model and the 
procedures to determine the existence of a QSH phase ap¬ 
plied to the homogeneous case, studying the relation be¬ 
tween interlayer coupling and the topological properties 
of the system. In section [HI] the same methodology is ap¬ 
plied to the case of a heterogeneous structure, graphene 
encapsulated by a trivial insulator, finding that topologi¬ 
cal properties can be induced even by trivial neighboring 
layers. Finally, in section IV we summarize our findings. 


II. HOMOGENEOUS MULTILAYERS 

Monolayer graphene consists of a triangular lattice 
with two atoms per unit cell that leads, in the recipro¬ 
cal space, to a hexagonal Brillouin zone that hosts Dirac 
cones in its corners. When N layers are considered the 
crystalline structure remains the same, only there will be 
2 N atoms per unit cell. We shall only use the so called 
Bernal stacking, shown in Fig. 0 , which is the ground 
state configuration, according to both DFT calculations 
and experimental evidence^-*®!. In Bernal stacked mate¬ 
rials an atom from the sublattice B(A) sits on top of an 
atom belonging to the other sublattice A(B). For N = 2 
there is only one way to achieve this, but for N > 2 there 
are different possible stacking orders. In figure [l] we show 
the different possibilities for N < 4, with a self-evident 
notation. 


A. The Model 

We describe the multilayers with the following tight- 
binding Hamiltonian: 

H = HmL + TjHi n te r + A L • S (1) 

where Hml and H inter account for the intralayer and in¬ 
terlayer hoppings, respectively, and the last term is the 


a) 



FIG. 1. a) Crystal structure of bilayer graphene system with 
a highlighted unit cell. Different colors for each layer are used 
to distinguish the two layers. In the inset the first Brillouin 
Zone is depicted with the high symmetry points and the Time 
Reversal Invariant Momenta colored in red. In b) Side view 
of the unit cells for all the different stackings studied. For 
the stackings with inversion symmetry the inversion center 
is shown at the crossing point of the dashed lines. For both 
figures red and blue denote sublattice. 


intra-atomic SOC. Our tight-binding model is based on 
four atomic orbitals, s,p x ,p y and p z . Both the intralayer 
and interlayer hoppings are described within the Slater- 
Koster formalism^. The intralayer hopping parameters 
are taken from ReP. In order to study the effect of 
interlayer coupling, the interlayer terms are scaled by a 
dimensionless parameter 77 . When 77 = 1, the ratio be¬ 
tween interlayer and intralayer V pp7T in graphene is taken 
0.13. Unless otherwise stated, in all our calculations 
we have 77 = 1. Within this model, the dimension of the 
Hilbert space for the minimal unit cell of the crystal with 
N layers is4x2x2x7V = 16V (4 orbitals per atom, 2 
atoms per layer, plus the two possible spin orientations). 

Without SOC, this model reproduces the very well 
known band structure of graphene (N = 1 ) and mul¬ 
tilayer graphene N > 1, that portraits these systems 
as zero-gap semiconductors. Within this model, SOC is 
known _to open a gap in the monolayeJ^ as well as in the 
bilayei® 39 52 . In the case of the monolayer graphene the 
gap is known to be topological. Within this model, the 
computed value of the gap lA6peV when we take a real¬ 
istic value of the atomic spin orbit coupling, A = lOmeV. 
This gap is much smaller than the ones obtained with 
accurate density functional theory (DFT) calculations, 
in the range of 30/ieV®. The reason for the discrepancy 
turns out to be that the mayor contribution to the SOC 
gap at the Dirac poi nt com es from the coupling to the 
higher energy d bands 29 52 . The later is a simple conse¬ 
quence of the fact that SOC opens a gap in second order 
in the coupling in the Dirac points when projected over 
the p band. In comparison, SOC acts as first order when 
considering channels involving the d band. Nevertheless, 
interlayer hopping may open a first order spin flipping 
channel in the p manifold, becoming of the same order as 
the intrinsic spin conserving d-level contribution. These 
last processes would be the ones missing in the multilayer 
Kane-Mele model, and should be added for completeness. 
In our case, for the sake of simplicity, we will focus on the 
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spin flipping channel, and use a four orbital tight-binding 
model considering A as a free parameter. 

Future work shall focus on the effect of the d-levels in 
multilayer graphene, which will not be addressed here. 

The effect of SOC on the band structure of the multi¬ 
layers can be summarized in the following points: 

1. SOC opens up a gap for all the N stacked layers 
considered, reproducing the existing results^ 9 ' for 
the case of N = 2. Notice that in the case of ABA 
and ABCB stackings, the system remains gapless 
until a critical value of A. This peculiarity is related 
to the non uniform evolution of the SO splitting of 
the linear and non-linear bands as shown in figure 

m 

2. The scaling of the gap with A is very similar for 
monolayer and N = 2, 3,4 multilayers as it is shown 
in Figure [2j Therefore, it is expected that within 
this model, the gap opened by the intrinsic SOC 
might be as small in multilayers as in monolayers. 

3. The magnitude of the band-gap is insensitive to 
the interlayer coupling. This result is somewhat 
surprising, since together with atomic SOC the in¬ 
terlayer coupling opens a spin-flip channel, other¬ 
wise missing in the monolayer case. In particular, 
switching on the interlayer coupling does not close 
the SOC gap of the monolayer as shown in Fig¬ 
ure [2| As a consequence, the ground state of two 
decoupled (77 = 0) monolayers can be adiabatically 
connected to the ground state of the bilayer (77 = 1). 

The last observation leads to the following result: odd 
N stacked graphene will be quantum spin Hall Insulators, 
whereas even N will not. More precisely, for a system of 
N decoupled monolayers the Z 2 invariant is: 

Z 2 (N) = [Z 2 (1)] n (2) 

Since the gap opened by A remains unaffected when 
switching on the interlayer coupling 77, the value of Z 2 
for graphene-like multilayers is also given by equation 
([ 2 ]). In the following we verify equation ([2| using two 
different strategies. In the case of inversion symmetric 
structures, we compute the Z 2 invariant. In all cases, 
we compute the edge states and check whether they fill 
the gap, or else. Independently on how the topological 
character is obtained, eq. § holds in all the cases. 

B. Calculation of the Z 2 invariant. 

Using the method developed by Fu and Kane in 2007= 
for systems with inversion symmetry it is possible to de¬ 
termine easily its topological character (the Z 2 invariant) 
by calculating the parity of the occupied Bloch wave func¬ 
tions at the time reversal invariant momenta (TRIMs). 

N 

; 0-1)" = IR 

m= 1 i 


where is the parity eigenvalue of the 2m th occu¬ 
pied state at the TRIM Ti = {T, Mi, M 2 , M 3 }. Using 
this method the topological character of a system will 
be determined just by the quantity ( — 1)^, resulting that 
(— l) v = +1 means trivial topology and (—l)*' = — 1 
means non trivial topology. The calculation for the sys¬ 
tems with inversion symmetry yields the following re¬ 
sults: 
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This guarantees that A and ABC crystals are topolog¬ 
ical but the bilayers and tetralayers (with inversion sym¬ 
metry) are not. This method cannot be applied to sys¬ 
tems without inversion symmetry, which are addressed 
in the next section using a different approach. 

C. Edge states 

To confirm equation § even for systems without in¬ 
version symmetry we look for the presence of gapless edge 
states. We consider armchair-terminated semi-infinite 
crystals. Using translation invariance along the direction 
parallel to the edge, we block-diagonalize the Hamilto¬ 
nian of the semi-infinite 2D crystal in terms of a collec¬ 
tion of fc|| dependent semi-infinite ID Hamiltonians, as 
indicated in figure ©• The ID Hamiltonian describes 
unit cells with AN atoms, where N stands for the num¬ 
ber of graphene layers. The intra-cell terms are denoted 
by Ho(k\\) and the inter-cell hoppings by V(k\\). 

The surface Green function of this block tridiagonal 
semi-infinite matrix can be written as: 

G ed 9 e (E, fc„) = [E + ie - ff 0 (fc||) - £ fl (fc||) - (fc„)] _1 

( 4 ) 

where E j R(fc||) is a self-energy that accounts for the cou¬ 
pling to the semi-infinite crystal, Y,H(k\\) is the self¬ 
energy due to its interaction with the H atoms included 
to get rid of the dangling bonds and e a small analytic 
continuation. 

The self-energy E# can be calculated employing a re¬ 
cursive Green’s function method that leads to the follow¬ 
ing coupled equations 

12 R (E,k n ) = V R (k ll )g R (E,k ll )Vl(k ll ) 

9R(E,k ||) = [E-Hoik^-XxiE,^)]- 1 (5) 

The E#(/c||) is calculated just as an additional itera¬ 
tion to the self-consistent calculation with the appropri¬ 
ate value for the hoppings C — H. 

For a given fc|| we compute the density of states using 

p{E,ku)~Alm[G ed ° e {E,k\\)} 


(3) 


(6) 
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FIG. 2. In a) the band structure close to the K point is shown for all the possible stackings of multilayer graphene with 
N = 2,3,4. Only when A ^ 0 (red line) a gap is opened at the Dirac points. Note that for ABA and ABCB stackings there are 
linear bands when A = 0 that when the SOC is switched on cause a smaller gap than in the other cases. In b) the dependence 
of the gap with the SOC A is shown. The anomalous behavior for the ABA and ABCB stackings is just due to the linear 
bands mentioned before. 



FIG. 3. Scheme of the mapping between a semi-infinite crystal 
and a semi-infinite chain. The coupling between each linear 
chain (with k\\ well defined) is introduced by means of a self 
energy Sj?. 


Using a similar approach we can also obtain the bulk 
density of states calculating the bulk Green function by 
recursion. 

In figure [4] we show the density of states for both bulk 
and edge for all the stackings as a contour plot in the 
k\\,E plane. For each stacking the left panel shows the 
bulk density of states, which are gaped for all the stack¬ 
ings and the right panel shows the edge states. The cal¬ 
culations are done for a rather large value of A = 2eV. 
The first thing to notice is that, for such large values of 
A, all the structures have edge states. However, only in 
the case of odd TV, shown in the left column, the in-gap 
states are gapless. This is a necessary condition in order 


to have a QSHI. In contrast, all systems with even N 
have edge states with a gap. Thereby, they are definitely 
not in the QSH phase, validating equation Q . There¬ 
fore, we conclude that odd N graphene stacks are QSHI 
and even N are trivial insulators. In all cases, the gap 
opened by SOC is quite small. 


III. HETEROGENEOUS MULTILAYERS 


In the previous section we have seen that for homoge¬ 
neous multilayers the gap opened by SOC has the same 
magnitude than for the monolayer. Thereby, homoge¬ 
neous multilayers of graphene would not improve the 
prospects for observation of the QSH phase compared 
to the monolayer. We thus explore the case of a het¬ 
erogeneous mult ilay er. This is motivated in part by re¬ 
cent experiments^ 7 that seem to indicate an enhancement 
of the SOC interaction in graphene due to proximity to 
WS 2 , a trivial semiconductor with quite large SOC and 
no inversion symmetry. There has also been plenty of 
work studying the enhancement of SOC interaction in 
graphene due to proximity to heavy metald^. However, 
it would be much more interesting if graphene could be 
driven into a QSH phase by proximity to an insulator, so 
that the only conducting channels would be only at the 
edges of graphene. 

Density functional calculations shov^ that a topolog¬ 
ical band-gap opens in graphene on top of both WS 2 
and WSe 2 , two widely studied two dimensional transition 
metal dichalcogenides (TMD). The magnitude of this gap 
is in the range of a few meV, i.e., two or three orders of 
magnitude larger than the intrinsic SOC gap. 

Here we propose a toy model to understand the open- 
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FIG. 4. For each structure bulk and edge density of states (left and right panel respectively). Gapless edge states appear only 
when an odd number of layers is considered independently of the stacking used. 


ing of a non-trivial gap due to proximity to a trivial in¬ 
sulator with strong spin orbit coupling. For that matter, 
we take graphene encapsulated between two monolayers 
of a trivial semiconductor with strong SOC and broken 
inversion symmetry. Specifically, the structure of these 
adjacent monolayers is that of a BN-like crystal (see fig¬ 
ure [ 5 ja)). The choice of the stacking is such that, glob¬ 
ally, the structure has inversion symmetry. Otherwise, a 
trivial band gap would be opened by proximity—. 

The BN-like crystal is described with the same inter¬ 
atomic Slater-Koster parameters than graphene, but very 
different on-site parameters. In particular we assume a 
large SOC A and a staggered potential ±m that breaks 
inversion symmetry of the top and bottom layers. Since 
we are interested in the proximity effect, we turn off the 
atomic SOC of the graphene layer. As in the case of 
the homogeneous multilayers, the interlayer coupling is 
characterized by the dimensionless parameter 77. In this 
case we impose zero SOC for the graphene layer, in order 
to study the proximity effect. For 77 = 0 the bands of 
this system would be the superposition of those of the 
top and bottom insulators, with gap 2m, and the bands 
of graphene, whose Dirac cones would lie inside the gap. 
Broadly speaking, this picture remains the same as the 
interlayer coupling is turned on. Interestingly, a non¬ 


trivial gap A opens in the Dirac cones only when 77 7^ 0 
and A / 0 . We have verified that this gap satisfies the 
scaling 


in the limit of small A, 77 and m -1 . This results implies 
that graphene can borrow SOC from a neighbor trivial 
insulator layer via interlayer coupling. Using the method 
of the TRIM we have verified that this insulator has Z2 = 
(— l) u = —1, and is therefore topologically non-trivial. 

The magnitude of the proximity effect away from the 
weak coupling limit of eq. 0 is shown in figures | 5 j We 
study the dependence of the proximity gap A as a func¬ 
tion of both the SOC A and the interlayer coupling 77 for 
two values of the encapsulating layer staggered potential 
m. It is apparent that, taking m = 2 . 0 eU (a trivial gap 
~ 1 . 5 eU) and A ~ 0 . 25 eV, values in line with those of 2D 
TMD, the proximity gap is in the order of ImeV, simi¬ 
lar to the DFT results. Therefore, our model provides a 
reasonable justification of the DFT computations, which 
are certainly more complete. 

Our toy model does not capture some probably im¬ 
portant features of real heterogeneous multilayers. For 
instance, the interlayer interaction could break inversion 
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FIG. 5. Panel a) shows the structure of the heterostructure 
considered. Panels (b,c,d) show the dependence of the in¬ 
duced gap in graphene due to the proximity of the encap¬ 
sulating layers. In panel b) it can be seen that the gap is 
proportional to A 2 and this estimation gets better as the gap 
of insulating layers gets bigger. Panel b) it can shows how 
the interlayer coupling 77 produces the expected effect, for 
small interlayer coupling the induced gap is small but it grows 
quickly as 77 increases. Panel c) shows the dependence of the 
induced gap with the sublattice imbalance 


symmetry which is expected to open a trivial gap. In 
addition, the geometry of our encapsulating layers was 
chosen to minimize the size of the unit cell, rather than 
to describe a real material. In general, the coupling of 


graphene to other 2D crystals will imply a new length 
scale, given by the size of the new unit cell. In this setup, 
the inversion symmetry breaking could average out. 


IV. CONCLUSIONS 


We have studied the quantum spin Hall phase in multi¬ 
layer graphene and in graphene encapsulated by a trivial 
semiconductor. In the case of multilayer graphene we 
find that only the stacks with an odd number of layers 
are quantum spin Hall insulators. However, the size of 
the gap is the same than for a monolayer, and thereby, 
most likely too small to be detected experimentally. In 
contrast, we propose a toy model for graphene encap¬ 
sulated between two semiconducting layers with strong 
SOC and a trivial gap. Our model shows that a non¬ 
trivial gap can be opened in graphene whose magnitude 
is controlled by the atomic spin orbit coupling of the 
adjacent layers. Our model provides a qualitative under¬ 
standing of recent DFT calculations 54 ^ as well as recent 
experimental work 37 and shows a promising route to ob¬ 
serve the quantum spin Hall phase in graphene. 
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